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ABSTRACT 

m ■ 

t— 1 \ We derive an analytic model for nonlinear "photon bubble" wave trains driven 

! \ by buoyancy forces in magnetized, radiation pressure-dominated atmospheres. 

>> \ Continuous, periodic wave solutions exist when radiative diffusion is slow com- 

; pared to the dynamical timescale of the atmosphere. We identify these waves 

with the saturation of a linear instability discovered by Arons — therefore, these 
wave trains should develop spontaneously. The buoyancy-driven waves are physi- 
cally distinct from photon bubbles in the presence of rapid diffusion, which evolve 
into trains of gas pressure-dominated shocks as they become nonlinear. 

Q-f Like the gas pressure-driven shock trains, buoyancy-driven photon bubbles 

can exhibit very large density contrasts, which greatly enhance the flow of ra- 
diation through the atmosphere. However, steady-state solutions for buoyancy- 
driven photon bubbles exist only when an extra source of radiation is added to 
the energy equation, in the form of a flux divergence. We argue that this term is 
required to compensate for the radiation flux lost via the bubbles, which increases 
with height. We speculate that an atmosphere subject to buoyancy- driven pho- 
ton bubbles, but lacking this compensating energy source, would lose pressure 
support and collapse on a timescale much shorter than the radiative diffusion 
time in the equivalent homogeneous atmosphere. 

Subject headings: accretion: accretion disks — hydrodynamics — instabilities — MHD- 
radiative transfer — X-rays: binaries 
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1. Introduction 

The tendency of "photon bubbles" to form in magnetized, radiation pressure-supported 
atmospheres is well-established. In photon bubble instability, the atmosphere spontaneously 
forms a propagating pattern of low-density channels separated by regions of high density. 
Radiation tends to diffuse through the underdense regions, avoiding the regions of high 
density. Both analytic calculations (Begelman 2001) and numerical simulations (Turner et 
al. 2005; Hsu, Arons & Klein 1997) indicate that the instability can lead to very large 
and persistent density contrasts. In this case, the atmosphere becomes "leaky" and the net 
radiation flux can greatly exceed the equilibrium value for a homogeneous atmosphere. 

Two families of linear photon bubble instabilities were discovered by Arons (1992) and 
Gammie (1998; see also Blaes & Socrates 2003), respectively. The Gammie modes depend on 
the interplay of gas pressure and radiation pressure and are destabilized slow magnetosonic 
waves. In the nonlinear limit they become trains of gas pressure-dominated shocks. The 
Arons modes are an outgrowth of internal entropy modes, and are related to previously- 
discovered instabilities driven by convection in a strong magnetic field, in the presence of 
heat conduction (Syrovatskii & Zhugzhda 1968). They are driven by buoyancy, are governed 
by the rate of radiative diffusion, and do not depend on the presence of gas pressure. The 
Arons modes occur where the parameter M , denoting the ratio of dynamical time to ra- 
diation diffusion time across the atmosphere, is smaller than one, indicating that radiation 
is relatively well-coupled to the gas. In contrast, the Gammie modes occur where radiation 
can diffuse across the atmosphere faster than the sound crossing time in the gas, a reflection 
of weak coupling — the relevant condition is M > /5 1//2 , where /5 <§C 1 is the ratio of gas 
pressure to radiation pressure. Both modes depend for their existence on the presence of a 
magnetic field stiff enough to enforce approximately one-dimensional motion. 

The purpose of this paper is to study the nonlinear development of the buoyancy-driven 
waves discovered by Arons (1992), and to assess their possible effects on radiation-dominated 
atmospheres. As in our earlier analysis of the nonlinear development of Gammie modes, we 
will use analytic techniques based on the assumption of periodic wave trains. In keeping 
with the essential physical characteristics of the Arons modes, we will ignore gas pressure in 
our analysis. This eliminates the possibility of obtaining a series of shock fronts, which prove 
to be the dominant feature of nonlinear Gammie waves. Instead, we will seek — and find - 
continuous periodic solutions. Because the waves are driven by buoyancy (and in contrast to 
our earlier analysis), we must carefully account for the vertical structure of the atmosphere. 
We do this by expanding our solution in powers of the vertical coordinate z, in addition 
to using a similarity variable to describe the waves. By focusing on wavelengths that are 
much shorter than the scale height of the atmosphere, we are able to perform a multi-scale 
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perturbation analysis in the radiation pressure (analogous to a WKB approximation, but 
nonlinear) to obtain a model for the nonlinear Arons waves. 

The plan of the paper is as follows. In § 2 we describe the basic equations and approxima- 
tions relevant to steady-state, nonlinear waves driven by buoyancy in a radiation-dominated 
atmosphere. The section concludes with a derivation of the basic equation describing the 
nonlinear waves. We solve this equation in § 3, where we also discuss the main features of the 
waves. In § 4 we interpret our results and discuss the waves' global effects on atmospheres. 
We emphasize that the existence of steady-state modes requires the presence of a distributed 
radiation source throughout the atmosphere. Absent this source, the atmosphere is expected 
to collapse due to the increasing leakage of radiation with height. We attribute the collapse 
of the atmospheres in the Hsu et al. (1997) simulations to this effect, and speculate that at- 
mospheres subject to buoyancy-driven photon bubble instability will generally lose pressure 
support and collapse on a timescale short compared to the radiation diffusion time in the 
equivalent homogeneous atmosphere. We summarize our results and present our conclusions 
in § 5. 



2. Equations 

We base our model on a simplified set of radiation-hydrodynamical equations, in which 
the magnetic field is stiff and vertical. The latter was assumed in Arons (1992), although 
Hsu et al. (1997) performed simulations for an inclined magnetic field. All results below can 
be readily generalized to an arbitrary (uniform) magnetic field direction. We assume that 
all vectors lie in the x — z plane, and that there is a uniform gravitational field —gz. 

The five equations that must be satisfied describe continuity, radiative diffusion in both 
the x and z directions, momentum conservation parallel to B, and conservation of radiation 
energy. Eliminating the radiation flux through the diffusion equation F = — (c//t)Vp/p, 
where k is the opacity (assumed constant) and p is the radiation pressure, we obtain (for a 
vertical magnetic field) 

p- t + (pv) z = 0, (1) 

V' t + vv z = -— - g, (2) 
P 



c 

3pt + 3vp~ z + Apv~ z = - 

K 



7/7 



- S , (3) 



where subscripts denote partial differentiation and dimensional coordinates are surmounted 
by a tilde. The source term —S'o in the energy equation represents an assumed divergence 
S'o = V • F of radiation flux that is intrinsic to the fluid. As we will show in § 3, this 
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added energy flux is necessary in order to obtain steady-state wave trains — otherwise the 
atmosphere collapses and no steady solution is possible. For simplicity, we will assume that 
So is a constant, corresponding to a linear variation of injected radiation flux with height. 
This is an arbitrary choice; steady wave train solutions may exist for various dependences of 
So on local state variables. 

We nondimensionalize the equations by normalizing the pressure, density, and velocity 
to fiducial values: p = pov, p = Pof]i v = v p u - We then define a pressure scale height 
H = Po/{pog) = c^/g, where Co is the isothermal sound speed associated with the radiation 
pressure. Dimensionless coordinates are given by x = x/H, z = z/H, t = v p t/H, and a 
characteristic Mach number is defined by m p = v p /cq. The diffusivity parameter from Arons 
(1992) is given by 

Mo = — (4) 
npoHco 

Finally, we define a dimensionless radiation source term by a = HSq/ '(4m p po c o)- The di- 
mensionless equations are then 

Vt + (v u )z = 0, (5) 
m 2 p (u t + uu z ) = - ^1 + ^ , (6) 

' -4a. (7) 



3z/ t + 3uu z + 4z/-u 2 = — - 



m p 



] + f — 

^ A v ^ 



We seek periodic "plane- wave" solutions with wavelength A, in which quantities associ- 
ated with the wave depend on position and time through the combination 

s = xt&n9 + z + t; (8) 

we henceforth denote differentiation with respect to s by a prime. Thus, the wavevector 
makes an angle 6 with respect to the z— axis, as in the analysis by Arons (1992). The 
quantity —v p is then the z— component of the phase velocity. The speed of the intersection 
of a wave front with the x— axis is given by — v p cot 9, while the overall phase speed of the 
waves is v p cos 9. 

In addition to the periodic dependence on s there must be a secular dependence on z. 
We will seek solutions in which the wavelength is short compared to the scale height, hence 
we can expand r/, v and u as a Taylor series in z, with \z\ <^ 1: 

z 2 

V = Vo(s) + zrj^s) + —r] 2 {s) + ... , (9) 



with similar relations for v and u. 
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Since we demand that all functions of s be periodic, we can define ^-independent wave- 
averaged quantities by the integral 



A(s)ds; (10) 

we are free to choose the density normalization so that (770) — 1- We seek steady-state 
solutions with no secular mass flow, implying that 



1 f nX 

lim — / r)u(ads + bdz) — 0, (11) 

n->-oo n\ ./n 



for all a and b, where 



z 2 

1]U = 1] U + z(7] 1 Uq + 7]qUi) + — (r) 2 U + 27^1 + 1] U 2 ) + . . . (12) 

This implies (r] u ) = and t]iUq + i] Ui = i] 2 Uo + 2rjiU\ + i] U2 = 0. Similarly, all higher 
order terms in the Taylor expansion for r\u must vanish. 

Let us now consider the continuity equation. To 0(z°) we have 

Vo + (VoUo)' + (Vi u o + V0U1) = Vo + (Vou )' = 0, (13) 
which is easily integrated to yield 

«„ = (14) 

Vo 

To 0(z l ) we have 

Vi + (Vi u o + V0U1Y + (V2U0 + 2r?iMi + r] u 2 ) =r][=0, (15) 
implying rji = —a = const., and 

ui = = a 5 — . (16) 

Vo Vo 

Note that the density scale height is H/a, with a — 1 corresponding to the p oc p oc 
atmosphere considered by Arons (1992). It is straightforward to carry the analysis of the 
continuity equation to higher orders in z with, for example, 772 = (3 = const, and u 2 = 
{2a 2 — /3t7o)(1 — Vo)/Vo- However, these higher-order terms are not required for our analysis. 



We next consider the energy equation. To O(z ) we have 



3(1 + Uq)v'q + 3-Uq^i + 4z/ (« + M i) 



M 



' M + "A + tan 2 ^ 



-4a, (17) 



- 6- 



while to 0(z l ) the equation is 

3(1 + Uq)v[ + 3mi(z/q + v-i) + 3m ^2 + 4z/ (-u' 1 + u 2 ) + 4z/ 1 (-u' + t^) 



m p 



V/i \VJ 2 \Vo Vo 



(18) 



where the terms in the expansion 



, v,/„ + HvA + TUj 2 + '" (19) 

are obtained from the momentum equation (6): 

= -1 - mj [(1 + UoK + «o«i] (20) 

= -Wp [(1 + uo)u[ + U U 2 + Ui(u' Q + Mi)] . (21) 

We argue below that {y z jr\) 2 is not needed at the desired level of approximation. Multiplying 
eq. (6) by rj, we also obtain the useful results: 

u 'o + v\ = -Vo ~ m l K + C 1 - Vo)ui] (22) 
v'x + "2 = - "^p K + ?7iMo + ^70^0^2] • (23) 

Our expansion in z is valid only if the wavelength projected on the z— axis is short 
compared to the scale height, A <C HcosO. (Note that this approximation just breaks 
down for the fastest growing linear waves, according to Arons [1992].) In this limit, the 
inertial terms of the momentum equation — those on the left-hand side of eq. (6) — should 
be small compared to the terms on the right-hand side. This ordering goes along with the 
assumption that the phase speed of the waves projected along the z— axis should be very 
subsonic, i.e., \m p \ <C 1. In this spirit, we introduce second expansion parameter for 

the dimensionless pressure, so that 

v = i/ + zvi + . . . = (1/00 + rripUQi + m 2 p u 02 + ...) + z(u 10 + m p u u + ...) + ... (24) 

We are free to set u 0Q = 1, v w = —1, in order to mimic an exponential atmosphere with 
p = poexp(—z/H), to lowest order. By renormalizing s to 

w = m^s (25) 

(where a prime henceforth denotes differentiation with respect to w), and treating all orders 
of u and r\ and their derivatives with respect to w as 0(1), we can expand equations (17) 
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and (18) in powers of m p . In these expansions we also treat M /m p and a as ~ 0(1). The 
leading terms in both equations are then ~ 0(m~ l ); in order to describe periodic wave trains 
we will need the 0(m~ l ) and 0(1) solutions to eq. (17), and the 0(m~ 1 ) solution to eq. (18). 

Carrying out this program, we first note that the terms proportional to (^ 2 /^)i in 
eq. (17) and (^2/^)2 in eq. (18) are ~ 0(m p ), so they can be ignored. Next, by expanding 
equations (22) and (23) we obtain: 

^01 = 1 - Vo] v' 02 = -(^11 + M o); "'a = ~(Vi + »2o)- (26) 

These results allow us to write down eq. (17) to 0{m~ l ): 

4u> = **> tan 2 e( U -k)' = ^ tan 2 6 ( = ^ tan 2 9 u> (27) 
m p \ 770 / m p \ VoJ m v 

(where we have used eq. [14]), implying 

m p = ^-Mo tan 2 9. (28) 

Now consider eq. (18) to 0(171^): 

4K-<) = ^tan 2 ^^-^V. (29) 
m p \Vo Vo J 

Substituting rji = —a and using equations (14), (16), and (28), we obtain v' u = r] — l = — 1/ 01 , 

where we have imposed the condition (i/ u ) = for a periodic solution. We may then set 

^11 = -^oi- 

Finally, we consider the 0(1) terms in eq. (17). Using previous results, we find that 
3(1 + Mo^oi + 3^0^10 — 0. Therefore, the surviving terms are 



4^01^0 + 4^1 = — 
m„ 



0°) +t an 2 0( — 
VoJ \Vo 



m p 



Mi ' ' x ' 



^ +tan 2 4^1^ 
Vo) V Vo 



-4a, 
(30) 

where we have used eq. (26) and substituted vu = Using the result that ^ 01 /?7o — u o 

and expressing u and U\ in terms of i] using equations (14) and (16), we obtain 

= sin 9 7. a, (31) 

MS Vo 

which is the equation describing steady wave trains. To simplify our study of solutions, we 
renormalize the independent variable to w = wsm9 and set a = aj sin 2 9, thus eliminating 
the factor sin 2 9 from eq. (31). The basic equation for our study of wave solutions is therefore 

vy = (i -■?)(«-?) (32) 

Tj J V 
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where a prime now denotes differentiation with respect to w and we have dropped the 
subscript "0" from i] . 

3. Solutions 

We seek continuous periodic solutions of eq. (32), subject to the constraint that (rj) = 1. 
That such solutions might exist is evident from the fact that eq. (32) can be written in the 
form 

where y is a function of rj and y' = rfdy/dr). f(y) is also an implicit function of rj through 
y. This is, of course, just a generalized version of the Hamiltonian equation of energy 
conservation in a potential. In order for eq. (33) to have periodic solutions, f(r)) must have 
two real roots, f(r] + ) = f(rj-) = 0, with / < for < rj < r] + , and the integral 

P + dy dri 

Jr, ^(-2/) 1/2 ' 

which equals half the wavelength, must be finite. 

The requirement of an intrinsic energy source in the atmosphere, a ^ 0, is also apparent 
from the form of eq. (32). If we choose a — 1, for example, then the first term on the right- 
hand side of eq. (32) is > 0, precluding a periodic solution unless a > 0. In general, for a 
given a and wavelength, a is an eigenvalue of the problem. 

Equation (32) is most easily solved by choosing a maximum density, r] + , then integrating 
to find the wavelength, minimum density rj-, and a, subject to the mean density constraint. 
It is better to start integrating from the maximum density rather than the minimum density 
because, in solutions with large density contrasts, i] + is exponentially sensitive to r\_. 

Figure 1 shows the density profile of the waves for r/ + = 1.1, 2, and 11, with a — 1. 
Parameters for a wider range of solutions are shown in Table 1. Even for wave trains with 
very low amplitudes, i] + — 1 1 (Fig. la), the density pattern is far from sinusoidal. At 
low amplitudes, r] + — 1 < 1, the depression in the low-density portion of the wave is sharper 
and deeper than the enhancement in the dense region. This behavior switches markedly for 
i] + — 1 > 1, and for high density contrasts (Fig. lb) the dense regions are extremely sharp 
and narrow, while the low density regions are broad and flat. 

Figure 2 and Table 1 illustrate the variation of minimum density wavelength A (in 
units of w), and intrinsic flux divergence a with \og w (r] + — 1). The minimum density varies 



(34) 
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in the opposite sense to the maximum density, but at large density contrasts this variation 
is slower than logarithmic. The wavelength decreases with increasing density contrast, also 
slower than logarithmically, while the intrinsic radiation source a, required for a steady-state 
solution, increases slightly faster than logarithmically, although still very slowly. 

Interestingly, the inverse relationship between density contrast and wavelength is op- 
posite to that found in the nonlinear shock train solutions (Begelman 2001). Numerical 
simulations (Turner et al. 2005) indicate that the shock train solutions evolve toward longer 
wavelengths, higher density contrasts, and larger luminosities. If higher density contrasts 
and larger luminosities are also energetically favored in the case of buoyancy- driven photon 
bubbles, then we would expect evolution toward shorter wavelengths. Indeed, the fastest- 
growing Arons modes have long wavelengths, close to the pressure scale height (see § 4.2), 
whereas the fastest-growing Gammie modes are short. 

3.1. Dependence on Density Gradient 

The parameter a represents the underlying density gradient in the atmosphere, with 
a — 1 corresponding to the "isothermal" atmosphere p oc p assumed by Arons (1992). The 
radiation entropy gradient Vln(p/p 4 / 3 ) ~ 4a/3 — 1, so all atmospheres with a > 0.75 are 
convectively stable according to the Schwarzschild criterion. Note that, because M < 1 
for atmospheres subject to buoyancy-driven photon bubbles, it may be possible for large 
convective cells to transport radiation energy without excessive leakage. (Recall that Mq 
is the ratio of the dynamical time to the diffusion time across a scale height.) The models 
with a = 1 and a = 2 correspond to convectively stable atmospheres, while the model with 
a = 0.5 is strongly unstable to convection. 

From Fig. 2b, it is apparent that density contrasts are relatively insensitive to a. The 
one glaring exception to this statement is the fact that no periodic solutions exist for low 
density contrasts, in the case of a — 0.5. This is easily understood by considering the right- 
hand side of eq. (32). In order to have a periodic solution, the right-hand side must be 
negative at the upper turning point (r) + ) and positive at the lower turning point (rj-). For 
sufficiently weak density contrasts and a < 1, it is possible to have a < rj- < 1. In this case, 
the first term on the right-hand side is negative, implying that a must be negative in order 
to make the right-hand side positive. But at r] + the first term is positive implying that a 
must be positive, leading to a contradiction. 

The quantities a and A are more sensitive to the density gradient than is r]-. At fixed r) + , 
a increases and A decreases with a. Since larger a means a steeper density gradient, these 
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a — 1 


a = 0.5 


a = 2 


logi (^+ - 1) 


?7_ 


AH 


a 


V- 


ah 


a 


V- 


ah 


a 


-2 


0.979 


82.2 


0.000111 








0.990 


6.29 


0.000146 


-1 


0.839 


25.9 


0.00979 








0.909 


6.24 


0.0136 





0.478 


9.44 


0.395 


0.308 


17.2 


0.386 


0.547 


5.18 


0.722 


1 


0.270 


6.11 


2.94 


0.225 


10.1 


1.49 


0.287 


3.97 


6.16 


2 


0.186 


5.39 


8.30 


0.170 


8.44 


3.79 


0.193 


3.62 


17.6 


3 


0.143 


5.11 


16.2 


0.134 


7.79 


7.43 


0.146 


3.49 


34.1 


4 


0.116 


4.96 


26.55 


0.111 


7.44 


12.2 


0.118 


3.41 


55.5 



Table 1: Parameters characterizing periodic solutions of the wave train equation (32), for 
three values of the density gradient parameter, a. All solutions satisfy the mean density 
constraint (77) = I. The wavelength A is in units of the dimensionless variable w. Solutions 
do not exist for a = 0.5 and r] + = 1.01, 1.1 — see § 3.1 for an explanation. 




Fig. 1. — Numerical solutions of the wave train equation (32) with a — 1, subject to the 
mean density constraint. Parameters for the plotted solutions can be found in the second, 
third and fourth rows of Table 1. (a) Solution with weak density contrast, 77+ = 1.1. (b) 
Superposed solutions with moderate (r? + = 2; heavy line) and large (r? + = 11; light line) 
density contrasts. 
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trends might reflect the increasing effectiveness of buoyancy in driving radiation through 
the atmosphere via the low-density channels. Larger a means that more energy is being 
transported by diffusion through the photon bubbles, while a smaller wavelength could mean 
that a smaller fraction of a pressure scale is sufficient to get the required driving force. Note 
that this is opposite to the trend for ordinary convection, in which energy is assumed to be 
advected by fluid motions, and energy transport is expected to be larger for smaller a. In 
buoyancy- driven photon bubble instability, diffusion rather than advection transports the 
energy. 



4. Interpretation 

4.1. Physical Properties of Solutions 

The dimensional wavelength of solutions, A, is related to the dimensionless wavelength 
A (in units of w), by 

X = m p H cot OX. (35) 
Combining this equation with eq. (28), we obtain 



— = -Motanfl. 
H 4 



(36) 



From our numerical solutions (Table 1) we have determined that A/4 ~ 0(1) for the cases 
with high density contrast. In order for our analysis to be valid, however, it is not sufficient 
for the wavelength merely to be smaller than the pressure scale height. A stronger condition 
- required by our expansion in powers of z — is that the wavelength projected onto the 




1 

0.8 
0.6 
0.4 
0.2 





-2-10 1 2 
log {T] + -1) 



Fig. 2. — (a) Plots of a and A vs. \og w (r] + — 1), for solutions with a — 1. (b) Minimum 
density r\_ vs. log 10 (?] + — 1) for a = 2, 1, 0.5 (upper to lower curve). 
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z— axis, \ z = A/ cos#, must be smaller than H. Thus, a necessary condition for the validity 
of our analysis is that 

cos 2 9 , , 

Mo < — -p 37 

From eq. (28) this condition implies that m p < sin6>/4 < 1, thus ensuring the validity of our 
expansion in m p . 

In all cases a is positive, implying that the atmosphere itself must supply extra energy at 
all heights, in order to support a steady-state wave solution. In physical units, the intrinsic 
flux radiated by the atmosphere is 

Fq ~ S H = 4<t sin 2 9 m p poc (38) 

Since o > O(10) for the high contrast solutions, these atmospheres must radiate furi- 
ously. By definition, the energy loss rate due to diffusion in the equivalent homogeneous 
atmosphere is ~ M p c . The inhomogeneous atmospheres in our model are losing energy 
asm 2 6* tan 2 9 times faster — possibly leaking energy on less than a dynamical time. 

It is important to remember that our adoption of an intrinsic flux divergence So was ar- 
bitrary — we took this step only because a steady-state wave solution could not be obtained 
without it. What would happen if, as is more likely, such a powerful intrinsic source term is 
absent? Although this case is not covered explicitly by our calculations, the only sensible de- 
duction would seem to be that the atmosphere loses radiation pressure support and collapses 
on the leakage time or the dynamical time, whichever is longer. Indeed, since the flux used 
in our radiation hydrodynamical equations is defined relative to the fluid frame, one way 
for the atmosphere to generate a flux divergence is for the gas to sink in the gravitational 
potential. We precluded this possibility by demanding (pv ) = 0, for mathematical reasons, 
but this outcome cannot be excluded physically. 



4.2. Relationship to Arons Instability 

While we have not proven that our nonlinear wave trains represent the outgrowth of 
Arons's (1992) instability, there are a number of aspects which make this identification 
plausible. First, our wave trains are driven by buoyancy forces — we guaranteed this by 
setting the gas pressure to zero (as did Arons in his linear analysis) and treating inertial terms 
as higher order than buoyancy terms. Second, we exploited the same near-cancellation in 
the flux-divergence as Arons did. In his slight imbalance drove the instability by 

diffusing more energy into the lower parts of bubbles than is removed at the top. In our 
case, it supplied the crucial driving term ij]' /rf)' in the wave equation (32). Third, our 
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waves only appear in the O(z) terms of an expansion in height, further illustrating that they 
are driven by buoyancy effects. Fourth, the vertical component of the wave phase velocity 
is negative, in agreement with the linear instability (Arons 1992, equations [19] and [42]; in 
contrast to the group velocity of the linear waves, which is positive — see Arons eq. [22]). 
Finally, the wave speeds in our model scale oc Mo, indicating that the wave speed is set by 
the diffusion of radiation across the bubble, as in the linear instability. 

According to the linear analysis, the most rapidly growing modes have kH cos 6 = 0.66 
(Arons 1992, eq. [44]), where k is the wavenumber. Projected on the z— axis, the wavelength 
is about 10 times longer than the scale height, putting it outside the limit of validity for our 
analysis. The fastest growing linear mode that is (marginally) consistent with our analysis 
has \ z ~ H. If we assume that this mode determines the orientation of the nonlinear wave 
pattern, then the inequality in eq. (37) is replaced by an approximate equality. In the limit 
M C 1, we then have cos 9 ~ /Mo <^ 1, and the dominant density contours will tend to 
align with the magnetic field. Thus, our nonlinear analysis reproduces another key feature 
of the linear instability, as well as the numerical simulations by Hsu et al. (1997). 

Finally, we note that the only existing numerical simulations of buoyancy- driven photon 
bubbles (Hsu et al. 1997), while reproducing the linear growth rates adequately, are plagued 
by the rapid collapse of the simulated atmosphere as the waves saturate. We conjecture that 
this outcome is real and reflects the lack of a compensating intrinsic source of radiation in 
the simulations. The fastest growing modes grow on a dynamical timescale (Arons 1992, 
eq. [45]), which is similar to the collapse time. 



4.3. Minimum Magnetic Field Strength and Effects of Gas Pressure 

The radiation pressure perturbations associated with the waves is Ap ~ (\/H)po tan 9 ~ 
M pc)tan#. These pressure fluctuations will be largely compensated by fluctuations in the 
magnetic field strength, in order to maintain overall pressure balance in the x— direction. 
The change in the magnetic pressure pb is related to sideways compressional or expansional 
motion by Aps ~ (Ax/X)pb, where Ax is the sideways displacement caused by the passage 
of the wave. Now, in order for the magnetic field to enforce approximately one-dimensional 
motion, we require that the sideways displacement be much smaller than the parallel dis- 
placement of a typical fluid element, Ax <C A/ cos9. This implies 

Pb A „ A , , sin 2 9 , , 

^» -sin0 = -Mo -. 39 

Po H 4 cost' 

Equation (39) appears somewhat different from the analogous criterion presented by Arons 
(1992) in connection with the linear instability. It is beyond the scope of this paper to carry 
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out a fully magnetohydrodynamical analysis. 

The effects of gas pressure will be most important in the densest regions of the waves, 
where 77 ~ rj + . Rapid energy exchange between matter and radiation will enforce approxi- 
mately isothermal conditions across the wavefronts, implying that the gas pressure is pro- 
portional to the local density. If the mean gas pressure is (3 g po, where /3 S < 1, then gas 
pressure will affect the wave dynamics when P g r] + p approaches Ap, or equivalently, when 

^-/^tanfl. (40) 

In such cases, gas pressure may limit the maximum density but is unlikely to have much 
effect on the values of A or a. 



5. Discussion and Conclusions 

We have produced an analytic model for periodic, nonlinear waves driven by buoyancy in 
a radiation pressure-dominated atmosphere with a stiff magnetic field. To isolate the effects 
of buoyancy, we ignored gas pressure and performed a multi-scale expansion in powers of the 
vertical coordinate z and wave Mach number m p . All solutions for steady-state, periodic wave 
structures require the existence of a large radiation source term in the atmosphere, which we 
modeled as a uniform, positive flux divergence term in the energy equation. Physically, we 
believe that this source term compensates for the enhanced loss of radiation via the waves. 
Without it the atmosphere would lose pressure support and collapse, on a timescale much 
shorter than the diffusion time associated with the equivalent homogeneous atmosphere. 

Based on several similarities of generic properties and physical mechanisms, we identify 
our wave solutions with the nonlinear development of the photon bubble instability discov- 
ered by Arons (1992). Simulations of the linear growth and saturation of this instability 
(Hsu et al. 1997) were plagued by the rapid collapse of the simulated atmosphere shortly 
after the waves reached nonlinear amplitudes. We suggest that this collapse resulted from 
the rapid release of energy via the waves, as indicated by our analytic model. 

Buoyancy-driven photon bubbles occur primarily when the radiation diffusion timescale 
across a pressure scale height is longer than the dynamical time, corresponding to the pa- 
rameter M being smaller than one. (It may be possible for buoyancy-driven photon bubbles 
to develop in atmospheres with M > 1, provided that the wavefronts are sufficiently oblique 
to the magnetic field.) This condition is satisfied in a number of astrophysical systems, 
including the equatorial regions of radiation-dominated accretion disks (where M ~ «, the 
viscosity parameter), accretion columns onto neutron stars in X-ray binaries, massive stellar 
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envelopes well below the photosphere, and gamma-ray burst fireballs. Whether the linear 
instability develops and saturates in any of these systems is not clear, because the growth 
of the instability is rather slow. For example, in accretion disks subject to the magnetoro- 
tational instability (MRI), the magnetic field structure is expected to change on a timescale 
comparable to the dynamical (Keplerian) time. This is shorter than the growth time of all 
but (possibly) the fastest growing buoyancy-driven photon bubbles. Hence, it would not 
be surprising if buoyancy-driven photon bubbles were suppressed in accretion disks. (Gas 
pressure-driven photon bubbles, on the other hand, can grow faster than MRI, and are ex- 
pected to be important in accretion disks [Begelman 2002, 2005; Turner et al. 2005].) In 
accretion columns onto a neutron star polar cap, however, the magnetic field is anchored into 
the neutron star crust and its structure is fixed. We would expect buoyancy-driven photon 
bubbles to grow to nonlinear amplitudes under these circumstances, perhaps leading to the 
collapse of the accretion column. 

What would happen to an atmosphere that collapsed under the action of buoyancy- 
driven photon bubbles? If we treat the collapsing atmosphere as a sequence of hydrostatic 
models, then both the diffusion parameter, M (oc Cq 1 for a fixed column density), and the 
ratio of gas pressure to radiation pressure, f3 g) are expected to increase. If M became large, 
then the atmosphere would enter the strong-diffusion regime and gas pressure-driven shock 
trains would become the dominant form of photon bubble. Alternatively, f3 g could increase 
until gas pressure halted the collapse. A third possibility, which may be relevant in the case 
of neutron star accretion columns, is that the gravitational binding energy liberated by the 
accretion flow provides at least part of the compensating flux divergence. It is not clear, 
however, whether this would be adequate to prevent the collapse of the column. 
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Theory Program grant NAG5- 12035 and by the National Science Foundation under Grants 
AST-0307502 and PHY-990794. Part of this work was carried out at the Kavli Institute for 
Theoretical Physics at the University of California, Santa Barbara; I thank the members of 
KITP for their hospitality. 
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